Three-dimensional dynamics of vortex lattice formation in Bose-Einstein condensates 
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We numerically study the dynamics of vortex lattice formation in a rotating cigar-shaped Bose- 
Einstein condensate. The study is a three-dimensional simulation of the Gross-Pitaevskii equation 
with a phenomenological dissipation term. The simulations reveal previously unknown dynamical 
features of the vortex nucleation process, in which the condensate undergoes a strongly turbulent 
stage and the penetrating vortex lines vibrate rapidly. The vibrations arise from spontaneous 
excitation of Kelvin waves on the proto-vortices during a surface wave instability, caused by an 
inhomogeneity of the condensate density along the elongated axial direction. 
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Atomic-gas Bose-Einstein condensates (BECs) provide 
a versatile testing ground for superfluidity, particularly 
in systems with an externally driven rotation. Several 
experimental groups have observed the formation of a 
lattice of quantized vortices in such a rotating BEC 
IH [H Q ■ The experiment that largely motivates the 
present study is a direct observation of nonlinear dy- 
namical phenomena such as vortex nucleation and lat- 
tice formation of a rotating BEC in real time |f| . There 
are some theoretical attempts to understand the dynam- 
ical properties of vortex lattice formation in this sys- 
tem H 0, H H [H [H 0. Among these studies, 
the experimental observations are well reproduced by 
two-dimensional (2D) numerical simulations of the time- 
dependent Gross-Pitaevskii equation (GPE) with a phe- 
nomenological dissipation term 0) B B E3 • 

This paper presents the full three-dimensional (3D) 
dynamics of vortex lattice formation in a rotating BEC 
through the numerical simulations of the time-dependent 
GPE with a dissi pati on term. Except for a few prelimi- 
nary works 0, IT]lll2l an d studies of the static configura- 
tions of vortices [14], HH , this problem has been analyzed 
using only 2D simulations. In contrast to 2D simula- 
tions, a 3D simulation gives significant improvement for 
describing real 3D systems. For the present study, these 
improvements are (i) the 3D simulation includes the axial 
(z) component, which increases the number of degrees of 
freedom of the fluctuation modes, and (ii) the 3D simu- 
lation can treat observed physical phenomena associated 
with "vortex lines", such as bending 14 fl, K elvin waves 
(helical displacements of the vortex core) | 16| , and recon- 
nection [ItJ ■ These features make the dynamical process 
of vortex lattice formation richer than that of the 2D case. 
Recently, Lobo et al. studied the 3D dynamics of rotat- 
ing atomic-gas BECs using classical field theory in which 
the dynamics is described by the energy-conserving time- 
dependent GPE and both the condensate wave function 
and its fluctuation modes are approximated by a sin- 
gle c- number field |11|. The dissipative process and the 
3D vortex dynamics of a rotating BEC in a spherically- 



symmetric trap was discussed in Ref. 12] . In our present 
study, we take account of the phenomenological dissipa- 
tion so as to reproduce the entire vortex formation pro- 
cess that was observed in the experiments [H 0] . 

In our previous numerical study E|, we described 
four dynamical stages of vortex lattice formation: (i) The 
condensate undergoes a large-amplitude quadrupole os- 
cillation that is damped because of the dissipation, (ii) 
After a few hundred milliseconds, the boundary surface 
of the condensate becomes unstable, generating surface 
ripples that propagate along the surface, (iii) The sur- 
face ripples develop into the vortex cores, which enter the 
condensate, (iv) The vortices form an ordered lattice and 
the global shape of the condensate recovers the axisym- 
metry. Particularly, in stages (ii) and (iii), the surface 
ripples are connected with proto-vortices, named "ghost 
vortices" 0, that appear in the low density region out- 
side the condensate. These processes are consistent with 
the observation in Ref. |5|, but the dynamics of the nu- 
cleating vortex lines themselves has not been reported 
yet. However, the above simulation was 2D, thus differ- 
ing from the actual 3D system. The present 3D simu- 
lations add some new qualitative aspects to the vortex 
dynamics, especially when the condensate undergoes a 
surface wave instability. 

Here we assume a cigar-shaped condensate because this 
is the case that was studied experimentally by Madi- 
son et al. 0| • In this study, we find that the tran- 
sitional process from the nonvortex state to the vortex 
state is governed by the strongly turbulent dynamics of 
the condensate wave function. This turbulence arises 
from spontaneous Kelvin wave excitation of the ghost vor- 
tex lines before the lines penetrate into the condensate. 
The Kelvin waves are excited because of the density in- 
homogeneity along the z-axis. Therefore, vortex lines 
also vibrate strongly immediately after they penetrate 
into the condensate. This significantly surpresses the 
visibility of the vortex cores in the 2D column density 
p(x,y) — J dzp(x,y, z) which is often used to visualize 
the vortices experimentally. 
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A BEC trapped in an external potential V(r) is de- 
scribed by a condensate wave function \&(r, t) obeying 
the GPE. In a frame rotating with a frequency f2 around 
the z-axis, the generalized GPE with a phenomenological 
dissipation constant 7 reads 0, H, M, LU| 

a* 



(a) 0.3 



-v 2 + ^ - M + u|*r (i) 



where /i is the chemical potential, L z — —i(xd y — yd x ) 
the angular momentum operator and the wave function is 
normalized to unity: J dr\^f\ 2 = 1. The units of energy, 
length and time are given by the corresponding scales of 
the radial harmonic potential as htuj_, b± = \Jh/2mu)x_ 



and 



respectively, where m is the atomic mass and 



u± the radial trap frequency. The trapping potential is 
V = [(1 + e)x 2 + (1 - e)y 2 ]/4 + A 2 z 2 /4 wi * n the aspect 
ratio A = uj z /uj± and a small anisotropy e; this form 
describes approximately the rotating potential used in 
the Madison et al. experiments PHH- The strength of the 
mean-field interaction between atoms is characterized by 
u, which is proportional to the s-wave scattering length a 
as u = 8irNa/b±, where N is the total particle number. 

We will compare the present results with 2D simula- 
tion results to clarify their differences. To reduce the 
dimension from 3D to 2D, we assume translation sym- 
metry along the z direction. Since we will be comparing 
to experiments on a condensate in an elongated, cigar- 
shaped trap with A ~ 10" 1 ~ 10" 2 0,0, we assume that 
the z-axis is along the symmetry axis of the experiments. 
We obtain the 2D formulation from the 3D approach by 
separating the degrees of freedom of the wave function 
as ^(r, t) = ip(x, y, t)<f>(z) [7|,|9j]. As a result, we obtain 
the 2D version of Eq. 0) with the alternating coupling 
constant u — > u 2 d = 87raA J dz\(f>{z)\ 4 / (J dz\<j)(z)\ 2 ) 2 ~ 
8iraNh z , where h z is a height of the cylinder. 

Because a vortex lattice state corresponds to a local 
minimum of the total energy functional, a phenomeno- 
logical dissipation constant 7 has been introduced in Eq. 
by replacing id/dt with (i - i)d/dt H Q3. We 
assume 7 = 0.03, which was obtained by Choi et al. 
by fitting a numerical simulation of the generalized GPE 
with the experimental data on collective excitations in 
Ref. [19j. Compared with the other terms in the GPE, 
this dissipation constant is so small that any small varia- 
tion of 7 could not change the dynamical process qualita- 
tively. Because the time development of the GPE with 7 
does not conserve the norm of the wave function, we ad- 
just the chemical potential fi at each time step to ensure 
normalization. This is done by calculating the correction 
A^ = (At)- 1 \n[f dr\t>(t)\ 2 / J dr\<S>(t + At)\ 2 } 22]. Since 
the dissipation term can be derived from the formulation 
for a finite-temperature BEC developed by Zaremba et 
al. |20j | (under some approximations ]^| ) , the dissipation 
should be related to the presence of a thermal component. 
Gardiner et al. obtained a similar but different GPE with 
the dissipative term using quantum kinetic theory |2l| . 
applying it to their simulation of vortex lattice formation 
[8j. For a fixed chemical potential, the dynamics in the 
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FIG. 1: (Color online) Time development of parameters in 
the 2D and 3D simulations, (a) The deformation parameter 
a. In both plots, the 3D results are the thick lines and the 2D 
results are thin lines, (b) The angular momentum per atom 
(L z ) (solid curve) and the total free energy F (dotted curve). 
The values of the parameters are A = 9.2 -1 , fl = 0.7u±, and 
U = 54760 for the 3D simulation and U2D = 440 for the 2D 
simulation. 



present study is exactly coincident with that in another 
model in Ref. Q provided that the atomic cloud and the 
trap are rotated at the same angular frequency and the 
gauge transformation $ — » i['e~ lAlt /'' is made. 

The numerical scheme to solve Eq. is a Crank- 
Nicholson method with the spatial mesh A x y z = 0.25 
and the time step At = 5 x 10~ 4 . We use the ground 
state solution of Eq. Q with a nonrotating axisymmetric 
trap (e = 0) as the initial state of the simulations. The 
conditions of the Madison et al. experiments on 87 Rb 
atoms are the following: JV = 3x 10 5 , ujj_ = 108.56 x 
2tt Hz, w 2 = 11.8 x 2vr Hz (A = 9.2" 1 ), and a = 5.29 nm, 
which yields the length scale &ho = 0.728 /an, the time 
scale ujJ 1 =1.47 msec, and the dimensionless coupling 
constant u = 54760. In the 2D calculation, we obtain 
the 2D coupling constant U2D = 440, where the height of 
the cylinder h z (along the z-axis) is assumed to be the 
Thomas- Fermi radius 2R Z = 180 /j,m A rotating drive 
is turned on by giving the value of fl = 0.7wj_ and non- 
adiabatically varying a small anisotropy parameter of the 
trapase(t) = AetQ(t g—t) with a rate Ae = 1.25X10 -3 / 
msec and a switch-off time t g = 20 msec, where 0(x) is 
a step function. 

To gain insight on the overall dynamics of vortex lat- 
tice formation, we first show in Fig. ^ the time devel- 
opment of the deformation parameter a = —Q((x 2 ) — 
(y 2 ))/(( x2 ) + (y 2 ))> the angular momentum per atom 
(L z ), and the total free energy F = (H — ilL z ) with 
H = -V 2 + V + m|*| 2 /2, where (A) = J dr**i*. 
From Fig. ^a), the deformation parameter initially has a 
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damped quadrupole oscillation as observed in the exper- 
iment 0- After two oscillation periods, a falls abruptly 
to a value just below 0.05 and settles down to a steady 
small-amplitude oscillation. Concurrently a vortex lat- 
tice forms during the rapid increase in the angular mo- 
mentum, and the simultaneous decrease in the free en- 
ergy as seen in Fig. ^b). This formation time, which 
depends on the value of 7, is in good agreement with the 
experimentally observed value of 250 msec [2j| . 

More detailed features are revealed in 3D images of the 
condensate density l^l 2 . The plots in Figure Eta) show 
the time development of a constant surface of the con- 
densate density. Before the abrupt decrease in a in Fig. 
da) , the condensate undergoes a quadrupole deformation 
[Fig. |2t a-i)]. After that, the condensate surface begins to 
show stripe-like surface waves that propagate along the 
surface and are nearly parallel to the z-axis [Fig. GJa-ii)]. 
The resultant surface waves lead to a state with turbu- 
lent density fluctuations as shown in Fig. ElViii). Most 
of these density fluctuations rapidly decay because of the 
dissipation, but some of them evolve into vortex lines 
that later penetrate into the condensate [Fig. EJa-iv)] 
and form a triangular lattice [Fig. EJVv)]. The surface 
wave instability occurs initially in the middle of the elon- 
gated axis (near z ~ 0). If the 3D condensate is regarded 
as an assembly of the sliced 2D condensates in the x- 
y plane, the central slice has the largest raduis and the 
largest 2D coupling constant i*2D • Then, if 2D analysis is 
applied, as in Refs. 0,13]' the plane with the largest cou- 
pling constant would be the first to have an instability, 
because the corresponding critical rotation frequency is 
smallest in that plane. This is consistent with the finding 
that the surface wave instability starts at z = 0. 

The most remarkable feature of the 3D simulations is 
that the vortex lines are bent. This feature can be seen in 
their column density profiles p{x, y) = J dz\^\ 2 shown in 
Fig. Efb), which have often used to visualize the vortices 
in experiments The column density profiles differ 

from the cross-sectional density at z = in Fig. Etc); in 
particular, the density dips characteristic of the vortex 
cores, the small dark circles in the image, are hard to 
discern just after the vortex penetration starts as shown 
in Fig. EJb-iii) and (c-iii). This is due to the bending 
of the vortex lines. After 350 msec, the blurred image 
of the vortex cores gradually becomes clear [Fig. Gfb- 
iv)]. We also found that the locations of the vortex cores 
in the cross-sectional density had an unsteady motion, a 
behavior that was not seen in the column density [25j. 
This motion indicates the existence of the Kelvin waves 
propagating along the vortex lines. Near equilibrium, the 
visibility of the vortices are still different from those in 
Fig. E£b-v) and (c-v) , because the vortex lines still bend 
near the edge of the condensates, which is consistent with 
the previous equilibrium arguments |14| . Animation of 
the above simulation is in Ref. plj . 

In the previous (2D) study 0, El , we investigated the 
dynamics of the phase field 8 = arg^. This is difficult 
in the 3D case; neverthless, the fact that the value of 



|v| diverges at the vortex core allows one to study the 
evolution of the phase field by calculating the superflow 
velocity |v| = |V0| = |(**V* - *V**)/2i|*| 2 |. Fig- 
ure Eta-i) and (b-i) shows that an array of vortex lines 
surrounds the condensate and precesses around it even 
before the surface wave instability occurs. These lines 
are ghost vortices that quickly appear near the ra- 
dial Thomas- Fermi boundary R±(z) after the rotation is 
turned on because their creation energy, which is propor- 
tional to the density, is almost negligible. The present 3D 
simulation reveals in Fig. EJa-ii) that the ghost vortices 
vibrate rapidly before they penetrate into the conden- 
sate. This is understood by the fact that the precession 
frequency of an off-centered 2D vortex is proportional to 
RJ 2 j2^j . By treating the 3D condensate as an assembly 
of 2D slices, one can understand that the precession of a 
ghost vortex near R±(z) at large \z\ is faster than that 
at z ~ 0. Therefore, a density inhomogeneity along the 
z-axis induces a mismatch of the precession frequency at 
each z, which causes the ghost vortex lines to twist and 
excites Kelvin waves. At the unstable stage, the vibrating 
ghost vortices give rise to violent density fluctuations in 
Fig. Eta-iii) , and the velocity field on the outskirts of the 
condensate enters the strongly turbulent regime. This 
feature appears to be unique to cigar-shaped conden- 
sates: spherically-symmetric condenesates and pancake- 
shaped condensates have nearly straight penetrating vor- 
tices After this turbulent stage, the dissipation re- 
moves some velocity fluctuations and makes the Kelvin 
modes damp out. At the quasi-stationary stage in Fig. 
Dfla-v) and (b-v) , the bent vortex lines inside the conden- 
sate move slowly and settle down to the lattice structure, 
while the surrounding ghost vortices precess and eventu- 
ally disappear outside. 

We briefly compare the 2D and 3D results. In Fig. 
da), the quadrupole oscillation frequency in 2D is the 
same as that in 3D, but the 2D result shows lingering 
damped oscillations after the surface wave instability. 
The 3D case is different because a given value of a in 
the x-y plane is integrated over z, and thus the oscilla- 
tion in a should vanish when the condensate density be- 
comes turbulent. Another difference between the 2D and 
3D cases is that the 3D system has a larger quasi-steady 
value of the angular momentum after about 0.3 seconds 
[Fig. njb)]. The reason for this difference is likely the 
following. After the vortex penetration, the centrifugal 
force expands the condensate in the radial direction. In 
the 3D calculation, the conservation of particle number 
decreases R z and results in an increase in the effective 
2D coupling constant U2D- In contrast, this coupling 
constant is fixed in the 2D case. The increase in U2D 
allows the penetration of more vortices , leading to a 
slight increase of the final angular momentum. There- 
fore, except for vibrations of the vortex lines, the overall 
dynamics observed in the experiments may be captured 
qualitatively in the 2D simulations. 

In conclusion, we did fully 3D numerical simulations 
of the GP equation to investigate the dynamics of vor- 



FIG. 2: (Color online) Time development of vortex formation in 3D. (a) The constant surface density |^| 2 = 2 x 10 -5 . (b) 
The column density p(x,y) — J dz\^(x,y, z)\ 2 . (c) The cross-sectional density \^(x, y, 0)| 2 . The times are (i) 95 msec, (ii) 220 
msec, (iii) 270 msec, (iv) 365 msec, (v) 500 msec. The box dimensions along x, y, and z axis are -14.5 to +14.5, -14.5 to +14.5, 
and -72.5 to +72.5, respectively, in units of feho = 0.728 [im (116 x 116 x 580 discretized space). 



tex lattice formation in a rotating BEC. We found that 
the generalized GP equation with a phenomenological 
dissipation term can accurately describe the experimen- 
tal observations reported in Ref. The 3D simula- 
tions revealed that, just before the vortex penetration, 
the condensate underwent a strongly turbulent process 
caused by the Kelvin wave excitation of the ghost vor- 
tices. These excitations were generated by a density in- 
homogeneity along the z-axis. The vortices were found 



to bend immediately after the ghost vortices developed 
into real vortices, which reduced the vortex visibility in 
the column density profile. 
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